I'm looking into doing a delta_sigma emulator. This is testing if the cat side works. Then i'll make an emulator for it.
In [1]:
from pearce.mocks import cat_dict
import numpy as np
from os import path
from astropy.io import fits
In [2]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
In [3]:
z_bins = np.array([0.15, 0.3, 0.45, 0.6, 0.75, 0.9])
zbin=1
In [4]:
zbc = (z_bins[1:]+z_bins[:-1])/2
print 1/(1+zbc)
In [5]:
a = 0.81120
z = 1.0/a - 1.0
Load up a snapshot at a redshift near the center of this bin.
In [6]:
print z
In [7]:
cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[a]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!
cat.load_catalog(a, particles=True, tol = 0.01, downsample_factor=1e-3)
In [8]:
cat.load_model(a, 'redMagic')
In [9]:
params = cat.model.param_dict.copy()
#params['mean_occupation_centrals_assembias_param1'] = 0.0
#params['mean_occupation_satellites_assembias_param1'] = 0.0
#my clustering redmagic best fit
params['logMmin'] = 12.386
params['sigma_logM'] = 0.4111
params['f_c'] = 0.292
params['alpha'] = 1.110
params['logM1'] = 13.777
params['logM0'] = 11.43433
print params
In [10]:
print cat
In [11]:
cat.populate(params)
In [12]:
nd_cat = cat.calc_analytic_nd(params)
print nd_cat
In [13]:
W = 0.00275848072207
print W
In [14]:
rp_bins = np.logspace(-1.1, 1.8, 11) #binning used in buzzard mocks
rpoints = (rp_bins[1:]+rp_bins[:-1])/2
theta_bins = np.logspace(np.log10(2.5), np.log10(250), 21)/60
tpoints = (theta_bins[1:]+theta_bins[:-1])/2
In [15]:
wt = cat.calc_wt(theta_bins, W, n_cores = 2)
In [16]:
import warnings
from scipy.interpolate import interp1d
rbins = np.logspace(-1.1, 1.8, 17) #make my own bins
xi = cat.calc_xi(rbins,do_jackknife=False,n_cores=2)
if np.any(xi<=0):
warnings.warn("Some values of xi are less than 0. Setting to a small nonzero value. This may have unexpected behavior, check your HOD")
xi[xi<=0] = 1e-3
#TODO I should be able to just specify my own rbins
rpoints = (rbins[:-1]+rbins[1:])/2.0
xi_rmin, xi_rmax = rpoints[0], rpoints[-1]
# make an interpolator for integrating
xi_interp = interp1d(np.log10(rpoints), np.log10(xi))
In [17]:
print np.log10(xi)
print xi_interp(np.log10(rpoints) )
In [18]:
import pyccl as ccl
# get the theotertical matter xi, for large scale estimates
names, vals = cat._get_cosmo_param_names_vals()
param_dict = { n:v for n,v in zip(names, vals)}
if 'Omega_c' not in param_dict:
param_dict['Omega_c'] = param_dict['Omega_m'] - param_dict['Omega_b']
del param_dict['Omega_m']
cosmo = ccl.Cosmology(**param_dict)
big_rbins = np.logspace(-1.1, 2.3, 41)
big_rpoints = (big_rbins[1:] + big_rbins[:-1])/2.0
big_xi_rmin = big_rpoints[0]
big_xi_rmax = big_rpoints[-1]
xi_mm = ccl.correlation_3d(cosmo, cat.a, big_rpoints)
xi_mm[xi_mm<0] = 1e-6
xi_mm_interp = interp1d(np.log10(big_rpoints), np.log10(xi_mm))
#correction factor
bias2 = np.power(10, xi_interp(1.2)-xi_mm_interp(1.2))
In [19]:
plt.plot(rpoints, xi)
plt.loglog()
Out[19]:
In [20]:
plt.plot(rpoints, xi)
plt.plot(big_rpoints, bias2*xi_mm)
plt.loglog()
#plt.xscale('log')
Out[20]:
In [21]:
plt.plot(rpoints, 10**xi_interp(np.log10(rpoints)))
plt.plot(big_rpoints, bias2*(10**xi_mm_interp(np.log10(big_rpoints))))
plt.loglog()
#plt.xscale('log')
Out[21]:
In [22]:
plt.plot(big_rpoints, xi_mm)
plt.loglog()
#plt.xscale('log')
Out[22]:
In [23]:
_theta_bins = np.radians(theta_bins)
tpoints = (_theta_bins[1:] + _theta_bins[:-1])/2.0
_wt = np.zeros_like(tpoints)
# need this distance for computation
x = cat.cosmology.comoving_distance(cat.z).to("Mpc").value#/self.h
In [24]:
assert tpoints[0]*x/cat.h >= xi_rmin #TODO explain this check
In [55]:
def integrand(u, x, t, bias2, xi_interp, xi_mm_interp):
r2 = u**2 + (x*t)**2
if r2 < xi_rmin**2:
return 0.0
elif xi_rmin**2 < r2 < xi_rmax**2:
return np.power(10, xi_interp(0.5*np.log10(r2)))
elif r2 < big_xi_rmax**2:
return bias2*np.power(10, xi_mm_interp(0.5*np.log10(r2)))
else:
return 0
In [56]:
from scipy.integrate import quad
large_scales, small_scales = [], []
for bin_no, t_med in enumerate(tpoints):
#log_u_ss_max = np.log(xi_rmax**2 - (t_med*x)**2)/2.0 #max we can integrate to on small scales
log_u_ls_max = np.log(big_xi_rmax**2 - (t_med*x)**2)/2.0 #max we can integrate to on small scales
'''
if not np.isnan(log_u_ss_max):
#print large_scales_integrand(log_u_ss_max, x, t_med, bias2, xi_mm_interp)
small_scales_contribution = quad(small_scales_integrand, -10, log_u_ss_max, args = (x, t_med, xi_interp))[0]
large_scales_contribution = quad(large_scales_integrand, log_u_ss_max, log_u_ls_max,\
args = (x, t_med, bias2, xi_mm_interp))[0]
else:
small_scales_contribution = 0.0
ls_start2 = big_xi_rmin**2 - (t_med*x)**2
log_ls_start = np.log(ls_start2)/2.0 if ls_start2 > 0 else
large_scales_contribution = quad(large_scales_integrand, log_ls_start, log_u_ls_max,\
args = (x, t_med, bias2, xi_mm_interp))[0]
'''
_wt[bin_no] = quad(integrand, -10, log_u_ls_max,\
args = (x, t_med, bias2, xi_interp, xi_mm_interp))[0]#(small_scales_contribution + large_scales_contribution)/cat.h #TODO check little h's
_wt*=W/cat.h
In [57]:
_wt
Out[57]:
In [58]:
np.radians(theta_bins)
Out[58]:
In [59]:
zbin = 1
wt_redmagic = np.loadtxt('/u/ki/swmclau2/Git/pearce/bin/mcmc/buzzard2_wt_%d%d.npy'%(zbin,zbin))
cov = np.loadtxt('/u/ki/swmclau2/Git/pearce/bin/mcmc/wt_11_cov.npy')
In [65]:
from astropy import units
def calc_wt(theta_bins, W, n_cores='all', halo = False, xi_kwargs = {}):
"""
TODO docs
"""
tpoints = np.radians((theta_bins[1:] + theta_bins[:-1])/2.0 )
wt = np.zeros_like(tpoints)
x = cat.cosmology.comoving_distance(cat.z)/cat.h
assert tpoints[0]*x.to("Mpc").value/cat.h >= rbins[0]
#ubins = np.linspace(10**-6, 10**4.0, 1001)
ubins = np.logspace(-6, 3.0, 501)
ubc = (ubins[1:]+ubins[:-1])/2.0
# TODO this is like this cause of a half-assed attempt at parraleization
# this is unesscary now, and could be discarded for a simpler for loop
'''
def integrate_xi(bin_no):#, w_theta, bin_no, ubc, ubins)
int_xi = 0
t_med = tpoints[bin_no]
for ubin_no, _u in enumerate(ubc):
_du = ubins[ubin_no+1]-ubins[ubin_no]
u = _u*units.Mpc/cat.h
du = _du*units.Mpc/cat.h
r = np.sqrt((u**2+(x*t_med)**2))#*cat.h#not sure about the h
if r > (units.Mpc)*cat.Lbox/10:
try:
int_xi+=du*bias2*(np.power(10, \
xi_mm_interp(np.log10(r.value))))
except ValueError:
#interpolation failed
int_xi+=du*0
else:
try:
int_xi+=du*(np.power(10, \
xi_interp(np.log10(r.value))))
except ValueError:
pass
wt[bin_no] = int_xi.to("Mpc").value/cat.h
'''
def integrate_xi(bin_no):
t_med = tpoints[bin_no]
int_xi = 0
t_med = tpoints[bin_no]
for ubin_no, _u in enumerate(ubc):
_du = ubins[ubin_no+1]-ubins[ubin_no]
du = _du*units.Mpc/cat.h
int_xi+=du.to("Mpc").value*integrand(np.log(_u/cat.h), x.value, t_med, bias2, xi_interp, xi_mm_interp)/cat.h
wt[bin_no] = int_xi
def integrate_xi2(bin_no):
t_med = tpoints[bin_no]
wt[bin_no] = quad(integrand, 1e-6, 1e3,\
args = (x.value, t_med, bias2, xi_interp, xi_mm_interp))[0]#/cat.h
#Currently this doesn't work cuz you can't pickle the integrate_xi function.
#I'll just ignore for now. This is why i'm making an emulator anyway
#p = Pool(n_cores)
map(integrate_xi2, range(tpoints.shape[0]))
#p.map(integrate_xi, range(tpoints.shape[0]))
#p.terminate()
return wt*W
In [66]:
#test integration methods
f = lambda x: x**2
def integrate_f():
xbins = np.linspace(0, 4.0, 501)
xbc = (xbins[1:]+xbins[:-1])/2.0
int_f = 0
for xbin_no, x in enumerate(xbc):
dx = xbins[xbin_no+1]-xbins[xbin_no]
int_f+=dx*f(x)
return int_f
def integrate_f2():
int_f = 0
return quad(f, 0, 4)[0]
print integrate_f(), integrate_f2()
In [67]:
wt2 = calc_wt(theta_bins, W, 2)
In [68]:
wt2
Out[68]:
In [69]:
tbins = (theta_bins[1:]+theta_bins[:-1])/2.0
#plt.plot(np.degrees(tpoints),_wt)
#plt.plot(np.degrees(tpoints), wt)
plt.plot(np.degrees(tpoints), wt2)
plt.errorbar(np.degrees(tpoints), wt_redmagic, yerr = np.sqrt(np.diag(cov)), fmt = 'o',
capsize = 50, label = 'Buzzard Mock', color = 'r')
#plt.plot(tbins, _wt)
plt.ylabel(r'$w_t(\theta)$')
plt.xlabel(r'$\theta$ [Arcmin]')
plt.xticks()
plt.loglog();
#plt.xlim([60, 240])
#plt.ylim([0, 2])
#plt.xscale('log')
Out[69]:
In [274]:
gt_data = hdulist[3].data
gt_rm, gt_bc = [],[]
for i, row in enumerate(gt_data):
if i == 20:
break
gt_rm.append(row[3])#gt_data[3,:20]
gt_bc.append(row[4])
In [ ]:
print gt_bc
print tbins*60
In [ ]:
plt.plot(gt_bc, gt_rm)
plt.plot(tbins*60, sigma_crit_inv*gt/cat.h)#/cat.h)
plt.ylabel(r'$\gamma_t(\theta)$')
plt.xlabel(r'$\theta$ [Arcmin]')
plt.loglog()
In [ ]: